Microscopic gauge-invariant theory of the c-axis infrared response of bilayer cuprate 
superconductors and the origin of the superconductivity induced absorption bands 
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We report on results of our theoretical study of the c-axis infrared conductivity of bilayer high-Tc 
cuprate superconductors using a microscopic model involving the bilayer-split (bonding and an- 
tibonding) bands. An emphasis is on the gauge-invariance of the theory, which turns out to be 
essential for the physical understanding of the electrodynamics of these compounds. The descrip- 
tion of the optical response involves local (intra-bilayer and inter-bilayer) current densities and local 
conductivities. The local conductivities are obtained using a microscopic theory, where the quasi- 
particles of the two bands are coupled to spin fluctuations. The coupling leads to superconductivity 
and is described at the level of generalized Eliashberg theory. Also addressed is the simpler case 
of quasiparticles coupled by a separable and nonretarded interaction. The gauge invariance of the 
theory is achieved by including a suitable class of vertex corrections. The resulting response of 
the model is studied in detail and an interpretation of two superconductivity-induced peaks in the 
experimental data of the real part of the c-axis conductivity is proposed. The peak around 400 cm~^ 
is attributed to a collective mode of the intra-bilayer regions, that is an analogue of the Bogolyubov- 
Anderson mode playing a crucial role in the theory of the longitudinal response of superconductors. 
For small values of the bilayer splitting, its nature is similar to that of the transverse plasmon of 
the phenomenological Josephson superlattice model. The peak around 1000 cm~^ is interpreted as a 
pair breaking-feature that is related to the electronic coupling through the spacing layers separating 
the bilayers. 

PACS numbers: 74.25.Gz, 74.72.-h 



I. INTRODUCTION 



The c-axis infrared response of the high-Tc cuprate 
superconductors (HTCS) is strongly sensitive to 
doping For underdoped HTCS, it reveals a surpris- 
ingly weak coupling between adjacent unit cells^ and a 
pronounced pseudogap.— In optimally doped materials, 
the real part of the normal state (NS) conductivity tJc 
is almost frequency- and temperature-independent for a 
broad range of frequencies and temperatures 4 In con- 
trast, the response of overdoped HTCS exhibits a metal- 
lic behavior."^ These findings, in particular the pseudo- 
gap, and the qualitative nature of the changes across 
the phase diagram, make the c-axis response one of the 
most interesting properties of the HTCS (for a review see 
Ref.i). In materials with two copper-oxygen planes per 
unit cell (the so called bilayer compounds), the c-axis re- 
sponse also reflects the electronic coupling within the pair 
of closely-spaced planes, that is of high interest for the 
following reasons: (z) Its renormalization with respect 
to the noninteracting case is an important fingerprint of 
the electronic correlations of the ground state, (ii) For 
underdoped HTCS, the manifestations of the pseudogap 
in (Tc interfere with those of the coupling. A prerequisite 
for an understanding of the c-axis pseudogap is thus a 
disentanglement of the former from the latter. (Hi) The 
couplingmay contribute to the condensation energy (see 
Refs. 0,11 and references therein). 



The character of the coupling has been debated since 
the early years of the high-Tc research. According 
to the conventional band theory, the hopping between 
the planes should lead to a splitting of the conduc- 
tion band into two branches: a bonding branch corre- 
sponding to states that are symmetric with respect to 
the mirror plane in the middle of the bilayer unit, and 
an antibonding branch corresponding to states that are 
antisymmetric.^ For some regions of the Brillouin zone 
(BZ), the bonding band is expected to be located below 
the Fermi level and the antibonding band above, which 
should give rise to the interband transitions.—^ 

The experimental normal state (NS) infrared spectra 
of the bilayer compounds, however, do not contain any 
structure that could be easily attributed to the transi- 
tions. Furthermore, the 20th-century photoemission ex- 
periments did not reveal the splitting of the conduction 
band. These findings could be interpreted in terms of 
strong electronic correlations localizing charged quasi- 
particles in individual planes, even in the case of the 
bilayer unit, and inhibiting the band splitting. The sim- 
ple band-structure based picture of the NS thus seemed 
to have failed. The experimental superconducting (SC) 
state infrared spectra of underdoped bilayer compounds 
exhibit features that are almost certainly related to the 
bilayer coupling: a broad absorption peak in the spectra 
of RecTc in the frequency region between 350 cm^^ and 
550 cm~^ (labeled as Pi in the following) and related 
anomalies of some infrared active phononSfii^iii These 
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features, however, also appear to be consistent with the 
absence of the conduction-band spHtting and the local- 
ization of charged quasiparticles: It was shown that they 
can be well understood and in some cases even fitted"'^^'^^ 
using the phenomenological model, where the stack of the 
copper-oxygen planes is represented by a superlattice of 
inter- and intra-bilayer Josephson junctions (the so called 
Josephson superlattice model, JSM).^'^ The mode Pi has 
been attributed to the transverse plasma mode of the 
model. A microscopic justification of the model in terms 
of quasiparticle Green's functions has been provided by 
Shah and MillisM. 

In the beginning of the 21st century, the situation 
changed. In particular, several groups have reported 
observations of two separate conduction bands in pho- 
toemission spectraJ^ii^iiiiiSiii^ The JSM is obviously not 
consistent with this observation. In addition, it became 
clear that the SC-state spectra of Reuc of YBa2Cu307_5 
(Y-123) exhibit two distinct superconductivity-induced 
modes: the mode Pi discussed above and another one 
around 1000 cm~^ (to be labeled as P2)M-^ It has been 
proposed, that the two are related, but in the light of the 
results of the recent systematic study by Yu et al.^'^ this 
appears to be unlikely. The presence of P2 cannot be 
accounted for in terms of the JSM. These facts thus call 
for a replacement of the simple phenomenological JSM 
with a more sophisticated theory involving the bilayer- 
split bands. Here we present such a theory and provide a 
fully microscopic interpretation of the superconductivity- 
induced modes Pi and P2. 

The basic ingredients of the theory are: (i) The lo- 
cal current densities, conductivities, fields, and a gener- 
alized multilayer formula. The local current densities of 
the intra- and inter-bilayer regions are expressed in terms 
of local conductivities and local fields. The fields differ 
from the average field because of charge fluctuations be- 
tween the planes. Macroscopic considerations of these 
charging effects lead to a formula for the total c-axis 
conductivity, that represents an extension of the com- 
mon multilayer formula^'^. (ii) The local conductivities 
are calculated using a microscopic model and the linear 
response theory. This is the main difference with respect 
to the phenomenological JSM, where they are estimated 
or obtained by fitting the data, (to) The microscopic 
description involves the two bilayer-split bands. The rel- 
evance of the bilayer splitting to the interpretation of the 
c-axis response has been pointed out in Ref. |2ll. (iv) The 
charged quasiparticles of the two bands are coupled to 
spin fluctuations. The coupling is treated at the level 
of generalized Eliashberg theory, as in Ref. [l^. (v) The 
gauge invariance of the theory, required for a consistent, 
i.e., charge conserving description of the charging effects, 
has been achieved by including a class of vertex correc- 
tions (VC) ensuring that the renormalized current ver- 
tices satisfy the appropriate Ward identities. The vertex 
corrections will be shown to lead to dramatic and quali- 
tative changes of the calculated response, similar to those 
occurring in case of the longitudinal response of a homo- 



geneous superconductor. 

Calculated spectra of Re CTc allow us to understand 
the nature of the peaks Pi and P2. The former will 
be shown to correspond to a collective mode resembling 
the Bogolyubov- Anderson mode of homogeneous super- 
conductors and the latter to a pair breaking (bonding- 
antibonding) peak. 

The rest of the paper is organized as follows. In Sec.HIl 
we present the essential aspects of the theory, the values 
of the input parameters and some computational details. 
Section Hill contains results and discussion. In Sec. IIII Al 
we focus on the relatively simple case of a BCS-like in- 
teraction between the quasiparticles. The analysis allows 
one to understand the consequences of the bilayer split- 
ting and the role of the vertex corrections, but the re- 
sulting spectra of Re CTc are not sufficiently realistic. The 
complex case of quasiparticles coupled to spin fluctua- 
tions is addressed in Sec. IIII Bl It will be shown that the 
calculated SC-state spectra display two distinct modes, 
similar to the experimental ones. Section fill CI presents 
a comprehensive discussion of the relation between the- 
ory and experiment including the interpretation of the 
superconductivity-induced modes. The summary and 
conclusions are given in Sec. IIVI The readers interested 
only in the main findings of the paper may consider skip- 
ping Sec. ini and some technical parts of Sec. IIII Al and 
SecHlBl 



II. THEORY 

In this section we elaborate on the basic ingredients 
of our theory mentioned in the introduction. First we 
briefly describe a phenomenological approach to the c- 
axis electrodynamics of the bilayer systems. In the sub- 
sequent paragraphs, we build up a corresponding micro- 
scopic description. 



A. Multilayer model 

The multilayer model proposed by van der Marel 
and TsvetkoAii^ provides a phenomenological descrip- 
tion of the c-axis electrodynamics of bilayer cuprates. 
These compounds are considered as consisting of homo- 
geneously charged copper-oxygen planes separated by in- 
trabilayer (bl) and interbilayer (int) spacing regions (see 
Fig. [1]). The dielectric function of the intrabilayer region 



eu{(^) =£00 + 



eouj 



(1) 



contains the interband dielectric constant Soo and the 
local conductivity au defined by jbi = o'bi^'bi, where jbi 
is the local current density and i?bi the local field. The 
interbilayer region is described in a similar way using 
the local conductivity tTint. To obtain the macroscopic 
(total) c-axis dielectric function £(w), modifications of 
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FIG. 1: (a) Crystal structure of Y-123. (b) Multilayer model, 
where intrabilayer and interbilayer current densities jhi and 
jint lead to a charge redistribution between the Cu02 planes, 
which modifies the local fields _Ebi and -Bint- 



the local fields due to the charging of the planes have to 
be considered. The result is 



d 



e{uj) Ehliuj) £int(w) 



(2) 



An extended version of the model, that we use in this 
paper, includes the dependence of the local current den- 
sities on both local fields: 



JL 



E 
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(3) 



The total c-axis conductivity crdw) is given as the ratio of 
the average current density (j) = (dbiibi + dintjint)/*^ to 
the average electric field {E) = (dbi-Ebi + dint^^int)/^- I3y 
employing the continuity relation between the charge and 
current densities jint — jh\ = dp/dt and the effect of the 
charged planes on the local fields, E^i — Eint = p/^o^oc 
(see Fig. [T]), we arrive at 



where 



f^blCTbl/bl + C^intfint/bl ^bl^bl/int + ^^intfint/int 



ibi + diata dua 



Eint _ CToo + Cbl/bl — CTint/bl 
Ehl CToo + CTint/int ^ O'bl/int 



(4) 



(5) 



and (Too = —iujeoSoc- The total dielectric function is 
given by e(tj) =£00 + iac{Lo)/eQUj. 

In the following we describe the calculations of the lo- 
cal conductivities aLL' based on a microscopic model. 
The subsequent incorporation of the interplane Coulomb 
interaction will then provide a microscopic justification 
for the phenomenological treatment of the plane-charging 
effects used in the model of van der Marel and Tsvetkov. 



B. Electronic structure — tight-binding bands, 
their renormalization and superconductivity 

One of the main components of our microscopic cal- 
culations are the two bilayer split bands. We therefore 
begin with the tight-binding description of these bands. 
The usual form of the in-plane dispersion 



-2i(cos kxa + cos kya) — At' cos kxa cos kya 



(6) 



will be considered, with the nearest neighbor and second 
nearest neighbor hopping matrix elements t and t' . The 
intrabilayer hopping is governed by the hopping matrix 
element tj^^, that is assumed to depend on and ky as 



t 



-Lfcii 



-(cos kxCi 



cos kyu)^ 



(7) 



This approximate form is suggested by the results of LDA 
calculationsS- and is roughly consistent with experimental 
data on Bi2Sr2CaCu208-H5 (Bi-2212)i^. Let us note that 
the essential results of our calculations do not depend 
on the form of i^fcy , what matters is the magnitude. In 
addition to the intrabilayer hopping, we consider a weak 
coupling through the interbilayer region with the matrix 



element of the same fc-dependence as ij 



The 



interlayer hopping splits the band 
bonding (B) and antibonding (A) - 



[6]) into two bands - 
with the dispersions 



2<±fe„^lfe„cos/c,d. (8) 



To account for the renormalization of charged quasi- 
particles and the superconducting pairing we adopt the 
spin-fermion model, where the quasiparticles are coupled 
to spin fluctuations. In the case of a single band, the 
model selfenergy (2 by 2 matrix) is given by 



j:{k,iE)^j^ J2 XSF{k-k',iE-iE')g{k',iE') 



k'AE' 



(9) 

which can be schematically written as the convolution 
E — g^XsF * G- Here g is the coupling constant, xsf is 
the Matsubara counterpart of the spin susceptibility and 
Q the Nambu propagator, Q{k,iE) — [iEtq — [ek — p)t^ — 
S(fc,ii?)]~^. The generalization to the two band case 
is straightforward and the selfenergies can be expressed 
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g Xsf ^va/b 



2 oven 

9 Xsf 



B/A 



(10) 



where we distinguish between the spin-susceptibility 
channels of even (xsF") odd (xsp*^) symmetry with 
respect to the mirror plane in the center of the bilayer 
unit. The diagrammatic representation of Sb/^ is shown 
in Fig. UKa). We have used the same form of xsf con- 
taining the resonance mode and a broad continuum as in 
Ref. [22 (details will be given in Sec. IIIF| ). 

Since the results of the selfconsistent calculations based 
on the spin-fermion model are difhcult to interpret, we 
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FIG. 2: (a) Diagrammatic representation of the selfenergies 
of the bonding (B) and antibonding (A) bands. The propaga- 
tors of the electronic quasiparticles and spin fluctuations are 
represented by the straight and the wiggly hues respectively, 
(b) Simple bubble approximation to the current-current cor- 
relator (|20p . Only the part given by Eq. (|2ip is shown. The 
black dots are the current vertices corresponding to j^^. (c) 
Current-current correlator with a renormalized current ver- 
tex (|23[) . The diagrams corresponding to the case of t'±k^^ ~ 
with no intraband contributions are shown, (d) Diagram- 
matic representation of the Bethe-Salpeter equation (|24|l . 
(e) Diagrammatic representation of the equation determin- 
ing the current-current correlator including plane-charging ef- 
fects. The dashed lines correspond to the interplane Coulomb 
interaction. 



first resort to the BCS level. The results obtained this 
way are easier to understand because of the absence of 
retardation and better possibiUties of analytical manipu- 
lations of the formulas. The even/odd interaction chan- 
nels are assumed to be equivalent which leads to the same 
superconducting gap = iAi„ax(cos fc^^a — cos kyu) in 
both bands determined by 



Afc = - 



E 

k\ne{A,B} 



A, 



kk' 



2E, 



■ tanh 



k' n 



(11) 



where V^fc' — —XwkWk' with Wk = {cos kxa — 
cosfcj,a)/2 is the BCS interaction of d-wave symme- 
try and E^A/B is the usual BCS quasiparticle energy 

For details see Ap- 



EkA/B 

pendix. 



C. Response to electromagnetic field 

Here we calculate the response of the model to the c- 
axis polarized electromagnetic field represented by the 
external vector potential A^^t = (0, 0, ^oxt) e"' ''^"*"'- 
The coupling of the tight-binding model to the elec- 
tromagnetic field can be obtained by multiplying each 
hopping term by the corresponding Peierls phase fac- 
tor according to the prescriptio n^^i^^i^^ '^'r'^r' ^ 

exp [-~{ie/h)Acxt ■ {R— R')] Cj^Cj^,. To fit the scheme 
of Sec. Ill A[ we formally distinguish between the vector 
potentials Au and A-mt, used for the hopping processes 
through the intrabilayer and interbilayer regions respec- 
tively. By expanding to the second order in the vector 
potentials, we arrive at the coupling Hamiltonian that 
can be used for extracting the c-axis paramagnetic and 
diamagnetic current density operators.^ The paramag- 
netic current density for q = 0, averaged over the corre- 
sponding region (bl/int), can be expressed as 



Jbl/int 



le 



Na?fi 



fc|l kzS 



Aks'^Bks 



"'^Bks^Akt 



'^k y-Aks^Aks ^Bks^Bks 



with the matrix elements 

2^±fc„(i±fc„ +t'±k„ cos k,d) 



J, 



(1) 



bl,fe 



iAk — ^Bk 



and 



(2) _ 2»^±fc|| ^Ifc 



L kzd 



£Ak — £Bk 



(12) 



(13) 



(14) 



The matrix element J^^^ is obtained from simply 
by interchanging tj_,^^^ and t'±k^^- 1^ *j_fc|| ~ ^ case, 

int.fc "^k 



where J^]''^ t^^ 



and jIVi. — j!^'' — 0, we arrive at 



the simplified expression 



ibl ~ 



Nna^h 



'^Akns'^Bk«s '^Sfcii s'^Afcii s 



fciis 



The summation runs over fc|| from the 2D Brillouin zone 
only and N is reduced accordingly. The diamagnetic cur- 
rent density is given by 



Jbl/int 



e^^bl/int ^bl/int 



(16) 



Cni,„ - C 



Aks^^Bk. 



Bks'-^Akt 



In the t'j^^ = case, Eq. (fTH]) simplifies to 



Jbl 



n-Bks — nAks 



(17) 



k„s 
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The total c-axis conductivity is constructed along 
the lines of Sec. Ill Al To this end, the current 
densities induced by the electric fields — + 
i6)AL {L E {bl, int}) have to be calculated and 
the local conductivities determined from jL{q,Lu) — 
'^LL'{q,i^)EL'{q,u;). At this point, the fields El' 
are still equal to the external field -Ecxt- However, it will 
be shown in Sec. Ill El that the local conductivities calcu- 
lated as outlined above, ignoring the charging effects play 
exactly the same role as in Eq. ([3|), i.e., they represent 
the response to the local fields. Within the framework 
of the linear response theory, the local conductivities are 
given by the Kubo formula 



in addition, an intraband contribution given by 



urj^\q 



0, ihv) 



k 



X Tv[gA{k,iE + iniy)gA{k,iE) 

+ gB{k,iE + ihv)gB{k,iE)] 



(22) 



This contribution has a similar frequency dependence 
as the in-plane conductivity, the main difference com- 
ing from the fe-dependence of the matrix element Jj^ 
Typically, it is rather small compared to (pij) . 



(2) 



{eyh^)KLSLL'+IlLL'iq,uj) 
i{uj -f iS) 



(18) 



D. Vertex corrections 



The first term in the numerator 

7(1) 



bl/int 



C^bl/int 



7Va2 



H Jh\/int,k^'r^Bks - riAks) , (19) 



k\\ kzS 



comes from the diamagnetic current densities and 
is related to the c-axis kinetic energy^^: In the 
= case, K\,\ = {d\,\/a?){T), where 
(T) is the intrabilayer kinetic energy per unit 
cell, T = -(l/iV||)X;fc||sij-fc||('^Bfc||s - nAk^^s) = 
Y.R^^ h;, . r;, (4r|| s^iR'^^ s + Cm,! .^2^. J • The 
second term in Eq. is the retarded correlation func- 
tion of the paramagnetic current densities 



^LL'iq,^) 



Na^dr' 



dte^-\[jl{q,t),jl,{-qM)0{t) 

(20) 

In the simplest approximation, the correlator is ob- 
tained by evaluating the bubble diagrams where the two 
current vertices are joined by two electron propagator 
lines. This is the approximation, where the vertex cor- 
rections are neglected. Since the propagators refer to 
the two bands, there are four possible combinations in 
total. Two of them correspond to interband transitions 
and their contribution to the Matsubara counterpart of 
pi)) at q = equals 



k.iE 



X T:r[gA{k,iE + ihv)gB{k,iE) 

+ geik, iE + ihiy)gA{k, iE)] (21) 

with the minus sign for L = L' and plus sign for L ^ L' . 
The corresponding diagrams are presented in Fig. [2fb). 
For t'_^f, 7^ 0, all the conductivity components acquire. 



The well-known deficiency of the simple bubble ap- 
proximations such as the one leading to Eqs. (PT|) and 
(|22| is the lack of the gauge invariance which mani- 
fests itself, e.g., by a violation of the normal-state re- 
stricted sum rule for the conductivity. For the normal 
state the conductivity components should satisfy the sum 
rule ^^ReaLL{oj)Auj = -iTie^ /2h^)KL. While the dis- 
crepancy between the left-hand side and the right-hand 
side in the corresponding case of the in-plane response is 
rather small (of the order of 1%^), here it is quite detri- 
mental - typically 20 — 30% - as demonstrated in Sec. lIIII 
Since there is an intimate relation between the gauge in- 
variance of the response functions and the charge con- 
servation, the large discrepancy indicates, that the conti- 
nuity equation between the current and charge densities 
is not even approximately satisfied. As a consequence, 
the use of the formula @ , which relies on the continuity 
equation, becomes questionable. In the following para- 
graph we show explicitly, how the requirement of gauge 
invariance enters a microscopic derivation of the formulas 
of Sec. 

To avoid the problems mentioned above, a gauge- 
invariant extension of the approximation (PT|) -)- (P^ is 
necessary. As found by Nambu^*, the gauge invariance 
of the response function is guaranteed if we replace the 
bare current-density vertex with a properly renormalized 
one. The required renormalization of this vertex (i.e., of 
the interaction of the quasiparticles with photons) is de- 
termined by the form of the quasiparticle selfenergy via 
the generalized Ward identity.^^ 

Here the situation is complicated by the presence of the 
two bands. To be able to express all the contributions 
in a systematical way, we first introduce the bare vertex 
factors {ie/Na^h)^^j^{k) (with m,n £ {A,B}) inferred 
from (|12p . In the corresponding diagram, the TO-th band 
propagator line with momentum k enters the current ver- 
tex of j£ and the n-th band propagator line leaves it. The 
possible combinations are: ~ ^7^a ^ '^^^^ 



_^int _ 7(1) bl 
lAB — •^int,fc' IBB 



•^bl.fc' IBA 



-lAA 



J, 



(2) 



The correlator Ull' involving the renormalized current 
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vertices r^,„(fc, iE, ifiv) 



fi^ Na?(i ^ 

^ k,iE,mne{A,B} 



Tr 



7^ 

' mn 



(k) 



xg^{k, iE)r^^{k, iE, ihi^)gn{k, iE + ihi^)] (23) 

contains two interband contributions with mn — AB and 
mn — BA. The corresponding diagrams are shown in 
Fig. mjc). For t'i^k^^ — 0, these are the only contributions. 

In the t', u ^0 case, also the intraband terms with mn — 
±fe|l r J 

AA, mn = BB contribute. 

The renormalized vertices T^^{k,iE ,ihv) consistent 
with the electronic selfenergies of the two bands obey 
the Bethe-Salpeter equations of the form diagrammati- 
cally shown in Fig.[5Jd). At this point, we have to distin- 
guish between the spin-fluctuation mediated interaction 
and the BCS interaction allowing for further analytical 
simplifications. Evaluating the diagrams in the former 
case we arrive at 

T'ABik^iE.ifiv) ^ -i^Asik) T0+ 
+ W ^ Xsr{k-k',^E-^E') 

^ k',iE' 

X (fc', iE')V\g {k', lE', inv)gA{k', lE' + ihv) 



f3N 



J2 xlf{k^k\iE^iE') 



k'AE' 



X GAik', z£;')r^^(fc', iE', inv)gB{k\ iE' + ihv) 

(24) 

and similar equations for the other renormalized vertices. 
Intraband current vertices and turn out to be 
simply the bare ones because of the symmetry of 7^^(fc), 
7^^(fc) (odd functions of k^) and g^-independence of 
XgF(q,a;) assumed in the t^^^^^ ^ case. In the BCS 
case, the interaction is non-retarded and separable, which 
leads to a simple fe-dependence and i_B-independence of 
F: r^g{k,ihv) = ^Asik) + XwkC^ {ihv). Here A is the 
BCS coupling constant and Wk is the d-wave symme- 
try function introduced in Sec. Ill Bl The Bethe-Salpeter 
equations and the current-current correlators can then 
be treated to a large extent analytically^, as shown in 
Appendix. In addition, the intraband contributions arc 
exactly zero in the optical limit of q ^ 0. 



E. RPA approximation of plane-charging effects 



In paragraph III Al we presented the results of a phe- 
nomenological approach to the effects due to the charg- 
ing of the planes. Here we outline a rigorous microscopic 
derivation of Eq. ([2]), where these effects are treated at 
the level of the random-phase approximation (RPA) . For 
the sake of simplicity, we restrict ourselves to the case of 
insulating interbilayer regions. 



The current density within a bilayer unit leads to a re- 
distribution of charge among the Cu02 planes. The elec- 
trostatic interaction of the corresponding charge densities 
is given by the interaction Hamiltonian 



Coulomb 



2eoo£o 



PP , 



(25) 



where p is the excess planar charge density. The current- 
current correlator Hbi/^i modified by this interaction 
along the lines of the RPA approximation, correspond- 
ing to the diagrammatic series shown in Fig. ^e) , reads 



H 



RPA 
bl/bl 



1 



eooSo + UP-p 



AJP- 



(26) 



where = Hti/bi, B^"'', B^'-J, and W'P are the 

correlation functions obtained without considering the 
charging effects. To proceed further towards Eq. 
we have to express these correlation functions using 
the conductivity-related current-current correlator H-'"-' 
only, eliminating H-'"'', IIP~^ , and H''"''. This can be 
achieved using the continuity equation for the charge and 
current densities. Let us note, that the gauge invariance 
of the local response functions is the necessary condition 
for the continuity equation to be valid. The result of the 
elimination can be written as 



^RPA 



O'bl/bl 



^0"bl/bl 



(27) 



The last step is the incorporation of the macroscopic 
averaging to obtain the macroscopic c-axis dielectric 
function 



eiu) 



equj (E) 



(28) 



where the symbols (j) and (E) denote the unit-cell aver- 
ages of the current density and the electric field, respec- 
tively. The averaged current density is given by (j) = 
{dhi/d)jhi, since the interbilayer regions are supposed not 
to contribute. The macroscopic field (E) consists of the 
homogeneous external field and the averaged field of the 
induced charge density {E) = i?oxt + (dhi/d) p/£o£oo. Us- 
ing the relation jbi — ''"S/bi-^cxt and the continuity equa- 
tion iujp = jbi, we obtain 



/■\ _ "^bl RPA p 
\J) - -^'^bl/bl^oxt , 



(E) = Se: 



d £00^0^ 



RPA 



E,. 



(29) 

Finally, by inserting these results in Eq. (|28p . we arrive 
at Eq. (121) with eint = £00 and Shi = £00 + iau/u/ so^- 
The local response function Cbi/bi calculated in Sec. Ill CI 
and III DI plavs the role of au- In the more general case of 
Eq. (|4]) , the derivation is analogous to the one presented 
here. We stress, that the use of Eq. ([21) or ([1]) is now 
accompanied by the requirement of the gauge invariance 
of the local conductivities. 
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F. Input parameters and computational details 

The values of most of the input parameters are the 
same as in Ref. [22- For the description of the bands 
we use the in-plane dispersion with t = 350 meV, t' = 
— 100 meV and the band filling n — 0.82. The values 
of the interplane hopping parameters will be specified 
later at the corresponding places in the text, since various 
regimes of the optical response corresponding to various 
values of these parameters are discussed. In the multi- 
layer formula, we use dbi — 3.4 A, dint = 12.0 A, i.e., the 
values corresponding to Bi-2212, and Soo = 5. 

The model spin susceptibility has the same form as in 
Refs. [13 and [22 containing the 40 meV resonance mode 
and a continuum with dimensionless spectral weights of 
0.015^./ and 0.016c respectively. In the t'^k^^ — ^ case, we 
distinguish between the channels of odd and even sym- 
metry and include the resonant mode with 6m = 1 in 
the odd channel only. The continuum with 6c = 2 is 
present in both channels. For t'j_^^^ ^ 0, the bonding 
and antibonding states are no more of the simple form 
\B),\A) = (|1)±|2))/V^, where |1) and |2) are state vec- 
tors residing on the first and the second plane of the bi- 
layer unit, respectively. The linear combination now con- 
tains fc-dependent coefficients. A proper construction of 
the the interaction vertices would extensively complicate 
the theory. To avoid this complexity, we take Bm = 1/2 
and 6c = 2 for both channels whenever t'^k^^ > ^■ 

The coupling constant .g = 3 eV was chosen to yield Tc 
around 90 K and the amplitude of the superconducting 
gap A around 30 meV. Some of the calculations were 
performed on the simpler BCS level, where we choose 
the value of the BCS coupling constant A leading to the 
same gap amplitude of 30 meV. 

The selfconsistent equations for the selfenergies ([TU]) 
and Bethe-Salpeter equations ((24)) were solved iteratively 
using a Brillouin zone grid of typically 64 x 64 x 32 points, 
and a cutoff of 8eV in Matsubara frequencies. In the case 
of small ij^max ^ 50 meV, the vertex corrections lead to 
a complete change of the response-function profiles and 
up to 10'^ iterations of the Bethe-Salpeter equation are 
required to achieve the convergence. The convolutions 
were performed using the FFT algorithm with the use of 
the symmetries of E and F. Since the calculations are 
very demanding in terms of computer time and memory, 
we have used gz-independent spin susceptibility which 
brings the advantage of fcx-independent E and F — 7. 
The calculated response functions were continued to the 
real axis using the method of Fade approximantsi^ 



III. RESULTS AND DISCUSSION 

A. Quasiparticles paired by the BCS interaction 

We begin with the simpler case of insulating spacing 
layers, i.e., t'j_ = 0. FigurelSfa) shows the local dielectric 
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FIG. 3: (a) Local dielectric function Ebi/bi in the simplest BCS 
case with t'±k^^ ~ 0, txmax = 45 meV, and Amax = 30 meV, 
vertex corrections are not included (NV). The thin (thick) 
lines correspond to the normal (superconducting) state, T = 
lOOK (T = 20K). The solid (dashed) lines represent the imag- 
inary (real) part, (b) The real part of the corresponding total 
c-axis conductivity obtained using Eq. ([2]). The thin (thick) 
line corresponds to the normal (superconducting) state. (c,d) 
The same as in (a,b) but with the VC included (for the su- 
perconducting state only). 



function £bi/bi of the intrabilayer region obtained using 
the bubble diagram of Fig.[21[b), i.e., with the vertex cor- 
rections neglected (this is abbreviated as NV). The thin 
(thick) lines correspond to the normal (superconducting) 
state, the solid (dashed) lines represent the imaginary 
(real) part. The NS response exhibits a sharp absorp- 
tion band near 80 meV due to the interband (bonding- 
antibonding) transitions. The SC-state response involves 
the superconducting condensate, which manifests itself in 
the real part of Ebi/bi and a pair breaking peak at llOmeV 
corresponding to final states with one Bogolyubov quasi- 
particle in the bonding band and one in the antibonding. 

Figure [3{b) shows the real part (Tc of the c-axis con- 
ductivity obtained using the multilayer formula [2] with 
£hi ~ £bi/bi and eint = £00. The dominant sharp peaks 
are located close to the frequencies of the zero crossings 
of 

£bi/bi- This can be understood using the fact that for 
Mbieintl < Mintewl Eq. ^ yields 



1 



Ctbltint 
dintEbl 



(30) 



and the expression on the right hand side has poles at the 
zero crossings of Recbi. Physically, the response is simi- 
lar to that of a system of thin metallic plates embedded in 
an insulating matrix, exhibiting a peak at the plasma fre- 
quency of the plates (the corresponding effective medium 
formulas can be found in Ref. |32| ). The narrow peak at 
160 meV of the NS spectra corresponds to the zero cross- 
ing of Re£bi associated with the interband transitions. 
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the peak of the SC state spectrum at 60 meV to the zero 
crossing due to the superconducting condensate. 

The VC change the response functions dramaticaUy: 
the SC-state spectrum of £bi/bi shown in Fig. ^c) dis- 
plays neither the superconducting condensate nor the 
pronounced pair breaking peak. They are replaced by 
a broad band centered at 70 meV. The real part of £bi/bi 
exhibits only two zero crossings (instead of the three oc- 
curring in the NV case, the difference is due to the ab- 
sence of the condensate). The one at lower energies is 
located in a region of strong absorption. As a conse- 
quence, the SC-state spectrum of CTc shown in Fig. [3]^d) 
displays only one pronounced maximum located at the 
same energy as that of the NS. 

Below we demonstrate, that the absence of the conden- 
sate in £bi/bi(VC) is a general consequence of the gauge 
invariance. The current density in the bilayer region in- 
duced by a homogeneous electric field E oriented along 
the c-axis can be expressed employing two different 
gauges of the electromagnetic potentials: 

(a) Aip = 0, Ec = iLoAc- Here Ec is the c-axis component 
of E^ is the scalar-potential difference between the 
two planes, and Ac is the c-axis component of the vector 
potential; 

(b) Ac = 0,Ec^ ~A^/du- 

Both approaches should lead to the same result. In 
the latter case the expression for the conductivity con- 
tains only a regular component proportional to a current- 
density correlator i^SiS The conductivity thus cannot pos- 
sess a singular component corresponding to the con- 
densate. Note that the above arguments utilizing the 
two gauges parallel those used when discussing the re- 
sponse of a homogeneous superconductor to a longitudi- 
nal probe. 

The analogy can be further used to understand the 
nature of the peak (mode) at 70 meV in Fig. [3jc) . We re- 
call that in homogeneous superconductors a longitudinal 
electromagnetic field excites the Bogolyubov-Anderson 
(B. A.) mode corresponding to density fluctuations of 
the electron system, associated with a modulation of the 
phase of the order parameter. "^"'■•^^ The energy of the 
B. A. mode is proportional to wf|q|, where vp is the 
Fermi velocity and q the wave-vector. So far we did 
not consider the Coulomb interaction between the carri- 
ers, that will shift the mode towards higher frequencies. 
In a single-layer superconductor (one Cu02 plane per 
unit cell), a longitudinal electromagnetic field with E \\ c 
would induce a B. A. -like mode with energy proportional 
to the Fermi velocity along the c-axis vpz, vfz ~ t±. In 
the present case of the IR response of a bilayer super- 
conductor the situation is more complicated. The elec- 
tromagnetic wave is transverse with q _L c. Nevertheless, 
it induces a charge density that is modulated along the 
c-axis. The modulation is analogous to the one associ- 
ated with the B. A. mode of a single-layer superconduc- 
tor with qr II c, |q| = 7r/d (d is the interplane distance). 
This is illustrated in Fig. ID The analogy allows us to 
interpret the mode as an analogue of the B. A. mode. 
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FIG. 4: (a) Schematic representation of the current density-, 
density-, and phase- pattern associated with the Bogolyubov- 
Anderson mode of a single- layer superconductor with q || c, 
|q| — Tv/d. (b) The same for the collective mode of the bilayer 
system discussed in the text. 
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FIG. 5: Dependence of the real part of the local conduc- 
tivity crbi/bi(r = 20 K) on the intrabilayer hopping ampli- 
tude txniax calculated with the vertex corrections neglected 
(a) and included (b). In (a), the absorption peak stops at 
2Amax ~ 60 meV when decreasing t±jnax- In (b), the energy 
of the peak is proportional to txmax. 



This point of view can be substantiated by comparing 
the Eqs. ((23l) . ((24|) with those describing the longitudi- 
nal response of a single-layer superconductor. For E || c, 
q II c, |q| — tt/cI, and for the Born-Karman region con- 
taining only two planes (a rather artificial situation), the 
latter possess the same form as the former. Note that 
the long-wavelength in-plane modulation of the electro- 
magnetic wave has qualitatively no impact on the mode. 

Figure [5] shows the t^-dependence of the intrabilayer 
conductivity (Jbi/bi calculated with the VC neglected (a) 
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and with the VC included (b). The frequency of the peak 
in (a) is determined by E^a + ^'fes which approaches 2A 
for t± 0. The energy of the collective mode in (b), 
however, does not depend on A; instead it is proportional 
to t± . This is consistent with the proposed interpretation 
of the mode and analogous to the relation uj{B. A.) ^ 
VF\q\. 

A further insight into the origin of the collective mode 
can be obtained by using arguments inspired by Ander- 
son's work on gauge invariance and the Meissner effect 
Anderson explains the difference between transverse and 
longitudinal excitations in terms of the complete second- 
order phonon-mediated interaction between electrons. 
The impact of the relevant interaction terms on longi- 
tudinal and transverse excitations is shown to be funda- 
mentally different. In the longitudinal case, these terms 
lead to a restoration of the gauge invariance, to the ab- 
sence of the condensate contribution in the response func- 
tion, and to the presence of a mode at a finite frequency 
proportional to the magnitude of the wavevector. In the 
present case, this role is played by the interaction terms 
involving the products of the form 



(a) 



and 



(31) 

They do not belong to the reduced BCS Hamiltonian 
leading to Eq. pT|) . They have, however, a profound 
impact on the final states since they provide an attrac- 
tive interaction between "elementary excited states", i.e., 
the states created by operators cj[,^|C^.g| and c^yijC^sx' 
that appear in the expressions for the current density 
operators on the right hand side of Eq. (fT^ . The result- 
ing bound state, i.e., the mode behind the maximum in 
the spectra, can be thought of as equivalent to a Cooper 
bound state of a pair of electrons - one from the bonding 
band and the other from the antibonding band - super- 
imposed on the BCS ground state of the two bands. 

It has been shown that the B. A. mode can be asso- 
ciated with oscillations of the phase of the order param- 
eter. We have checked that for small values of t± the 
collective mode of our bilayer case can be similarly asso- 
ciated with oscillations of the relative phase of the two 
planes. The pattern of the phase modulation is shown 
in Fig. m Finally, the mechanism of the increase of the 
frequency of the mode when going from the local con- 
ductivity dbi/bi to the total conductivity, involving the 
Coulomb interaction of the charged planes, is an ana- 
logue of the Anderson-Higgs mechanism. 

Next we address the more complicated case of t'^ ^ 0, 
where the theory involves the four local conductivities 
defined by Eq. HH): (Tbi/bi, crbi/i„t, (Ti„t/bi (that differs 
from o-bi/int only by a factor of dbi/^int), and (Tint/int- 
Figure inj a) shows the real parts of Cbi/bi, fbi/int, and 
fint/int for representative values of the hopping parame- 
ters. The dashed (solid) lines correspond to the NV ap- 
proximation (to the approach with the VC included) . In 
the NV case, all the conductivities display a pronounced 
structure around 100 meV: a maximum in RecTbi/bi and 
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FIG. 6: (a) Real parts of the local conductivities (Tbi/bi, 
''"bi/int, and o"int/int Calculated considering the BCS interac- 
tion between the charged quasiparticles for t±max = 45 meV, 
4max = 0.2fxmax, A^ax = 30 HieV and T = 20 K. The spec- 
tra of Reubi/int and Reaint/int are four times magnified. The 
dashed (solid) lines correspond to the NV approximation (to 
the approach with the VC included), (b) Real part of the 
total c-axis conductivity calculated using Eq. Q for various 
values of the ratio fYmax/^Xmax. The figure demonstrates the 
effect of the interbilayer hopping. 



RetJint/inti and a wave- like feature in RetTbi/int- The VC 
lead to drastic changes of Re (Tbi/bi- The maximum shifts 
towards lower energies and its spectral weight increases 
on the account of the condensate (not shown). On the 
other hand, the structures in Re (Tint/int and Redbi/int re- 
main qualitatively the same and, in particular, they do 
not shift towards lower energies. 

The difference can be understood using Fig. [T] Part 
(a) provides a schematic representation of the Wannier- 
like orbitals of the planes and of the interplane hopping 
processes. For t'j_ <^ t±, it is useful to consider bonding 
and antibonding orbitals of the individual bilayers shown 
in (b), 4-5 = (1/V2)[«'i + 4-2], ^A = il/V2)[^i - 4-2]. 
The local current densities and conductivities can be dis- 
cussed and understood in terms of the transitions de- 
noted by the arrows. The intra-bilayer current-density 
operator is connected with transitions within individual 
bilayers, marked by the solid arrows. Note that these 
transitions create two quasiparticles from the same bi- 
layer unit. The inter-bilayer current-density operator 
is connected with transitions between adjacent bilayers, 
marked by the dashed arrows. These transitions create 
two quasiparticles from different units. The important 
point is that the final states due to the former (latter) 
transitions are strongly (weakly) modified by the VC be- 
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FIG. 7: Scheme illustrating the differences between the local 
conductivities (Tbi/bi and (Ji-nt/int- The Wannier-like orbitals 
of the planes are denoted by and ^2, the bonding and 
antibonding orbitals of the individual bilayers by ^'a and 4' b 
respectively. The current density operators j'^i and j?^^ can 
be associated with transitions marked by the solid and the 
dashed arrows respectively. 



cause the interactions are restricted to individual bilayers. 
This is the reason, why Cbi/bi (determined by the matrix 
element of Eq. (^0]) involving states generated by j^j, i.e., 
by the former transitions) is strongly influenced by the 
VC, whereas Uint/int (determined by the matrix element 
of Eq. ([20|) involving states generated by j-Jj^, i.e., by the 
latter transitions) hardly changes. The changes of cTbi/int 
are more complex because different final states are in- 
volved: one due to the former transitions and the other 
due to the latter transitions. 
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FIG. 8: Dependence of Reabi/bi (a) and Recrint/int (b) on 
the hopping amplitude f±max. A constant value of the ratio 
*Xmax/^Xmax of 0.2 has been used. If the vertex corrections 
are included, the B. A. -like mode characterized by a linear 
tx-dependence appears in the bl/bl component and, for small 
values of txmax, also in the int/int component. 



B. Quasiparticles coupled to spin fluctuations 



The t_L-dependencies of RetTbi/bi and RecTint/int are 
contrasted in Fig. [51 It can be seen that the VC cause 
qualitative changes of RcCTbi/bb in particular, they lead 
to the linear ij^-dependence of the frequency of the max- 
imum. The spectra of Reaint/intj on the other hand, do 
not change qualitatively, except for the very small values 
of i^, where the B. A. -like mode appears even in the 
intcrlaycr conductivity. 

The structures of the conductivities (Tint/int and crbi/int 
give rise - via the multilayer formula - to a maximum 
in the spectra of Re Cc as shown in part (b) of Fig. [6] 
It can be seen that the magnitude of the peak is pro- 
portional to t^. The maximum of Re (Tint/int and the 
related peak of Re tXc can be interpreted simply as an in- 
terband bonding-antibonding pair-breaking (coherence) 
peak, with the coherence factor proportional to the mag- 
nitude of the band splitting. We recall that in one-band 
superconductors, the conductivity does not exhibit any 
coherence peak around 2A due to the fact that electro- 
magnetic absorption belongs to phenomena governed by 
case II coherence factors Here the situation is differ- 
ent because of the presence of the two bands. 



Here we discuss the c-axis response obtained by 
the selfconsistent computations within the spin-fermion 
model. As in the previous section, we begin with the 
simpler case of t'^ = Q. Interestingly, the VC play an im- 
portant role even in the normal state, an effect that was 
not addressed at the BCS level. In particular, for the NS 
local conductivity calculated with the VC neglected, the 
restricted conductivity sum rule 

I,^— Re (Tbi(w) dec; = -i^bi (32) 
7i"e^ Jq+ 

is strongly (by tens of percent) violated. This is demon- 
strated in Fig. [9l It can be seen that the deficiency is 
more serious for smaller values of the band splitting. 
With the VC included, the sum rule is satisfied. Let 
us note that the corresponding violation of the in-plane 
sum-rule is an order of magnitude weaker, 

Changes of CTbi/bi caused by the incorporation of the 
VC are so pronounced that there is not much similarity 
between the NV and the VC spectra, see Fig. [TUl The 
NS spectra for the NV case shown in (a) display a broad 
low-energy absorption band corresponding to bonding- 
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FIG. 9: Optical weight h = Ih^ jixe^ /J^ Reabi(w) dw of the 
local conductivity (Jbi divided by \K\^\ \ as a function of temper- 
ature for various values of the interplane hopping parameter 
fx. Here K\a is the effective c-axis kinetic energy defined by 
Eq. (|19|l . The results were obtained using Eq. (|2ip within the 
spin-fermion model with no interbilayer hopping (txmax = 0) 
and the resonant mode in the odd interaction channel only, 
the vertex corrections have been neglected. With the inclu- 
sion of the vertex corrections, Jo /|Js'bi| is always equal to 1, 
both for the normal and the superconducting states. 

antibonding interband transitions, that can be compared 
with the sharp absorption structure of the NS spectra of 
the BCS case. The broadening with respect to the BCS 
case is due first to the finite lifetime of the quasiparticles 
and second to the presence of a pronounced incoherent 
background of the spectral function. The VC shift the 
absorption band towards higher energies, see part (b). 
This can be understood in terms of the complete second- 
order spin-fluctuation-mediated interaction between the 
electrons: the relevant terms can be shown to correspond 
to a repulsive coupling of the excited states with one 
electron in the antibonding band and one hole in the 
bonding band. In the total c-axis conductivity shown in 
Fig. fTUT e') the band is shifted even further and it is very 
broad. 

The spectra of the SC-state for the NV case as shown 
in Figlin] (c) display a pair breaking peak at about 2A, 
corresponding to a similar feature of the BCS case, and 
a continuum with an onset around 80 meV. In addition, 
Reubi/bi also contains the contribution of the condensate 
~ S{lo) (not shown). The VC transform the spectra in a 
similar way as in the nonretarded case, see part (d). They 
destroy the condensate and the pair breaking peak; in- 
stead, a sharp maximum (mode) appears, whose energy is 
proportional to As discussed in Sec. IIII Al the mode 
can be interpreted as an analogue of the B. A. mode. 
In the following we shall call it simply B. A. mode. In 
the total c-axis conductivity, see Fig. [TUTf). the mode is 
shifted towards higher energies by virtue of the Coulomb 
effects associated with the charging of the planes, the cor- 
responding peak will be labeled as T\. Let us emphasize 
that the sharp peak shows up only in the superconduct- 
ing state, the presence of a narrow mode in xsf is not a 
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FIG. 10; Effect of the vertex corrections on the conductiv- 
ities. Real part of the local conductivity cTbi/bi [(a')-(d)] and 
of the total c-axis conductivity [(e)-(f)] calculated within the 
spin-fermion model. The spacing layers have been taken to 
be insulating, i.e., t^max = and the neutron resonance has 
been included in the odd-interaction channel only. Results for 
several values of the hopping parameter txmax are presented. 
The spectra in (a) and (c) [(b), (d), (e), (f)] have been ob- 
tained for the NV [VC] case. The thin line in (f) represents 
the estimated energy dependence of the spectral weight of the 
peak labeled as T\. For txmax = 30 meV (45 meV, 60 meV), 
the estimated value of the spectral weight is 4000 f2~^cm~^ 
(4600n~^cm-^ 2300 n-^cm"^). 



sufficient condition for its appearance. 

In the remaining part of this subsection, we address 
the case of nonzero t'^. Similarly as in the BCS case, 
the impact of the VC on (Xbi/bi is much stronger than 
that on the other local conductivities. This is illustrated 
in part (a) of Fig. [TlJ The NV spectra of RecTbi/bi ex- 
hibit a pair-breaking peak at about 2 A, similar as in 
Fig. [TUTc). The VC destroy the peak and lead to the 
formation of the B. A. mode. For the present value of t_L 
of 150 meV, the mode is located in the region of the con- 
tinuum [cf. Fig. rTUT d)] and thus only gives rise to a weak 
structure around 100 meV. The inter-bilayer conductiv- 
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FIG. 11: Effect of finite interbilayer hopping (txmax > 0) 
on the conductivities. Real parts of the local conductivities 
o"bi/bi and cTint/int (a) and of the total c-axis conductivity (b) 
calculated within the spin-fermion model. The neutron res- 
onance has been distributed equally in both (even and odd) 
interaction channels. The dashed (solid) lines in (a) represent 
the spectra obtained with the VC neglected (included). The 
spectra in (b) have been obtained with the VC included, the 
dashed (solid) lines correspond to the normal (superconduct- 
ing) state, results for several values of txmax are shown. 
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FIG. 12: (a) Energies of the B. A. peak and of the 2A-peak 
in the local conductivity (Tbi/bii and of the structures Ti and 
T2 in the total c-axis conductivity (for examples, see Figs. 1101 
and lll|) as a function of the intra-bilayer hopping parameter 
txmax. A constant value of the ratio ^xmax/^J-max of 0.2 has 
been used, (b) Real part of the total conductivity for values 
of txmax, where both structures Ti and T2 are visible. As 
ixmax decreases, the peak Ti emerges from the background 
at txmax ~ 100 meV, and it quickly becomes the dominant 
feature. Eventually, it covers the r2-peak at txmax ~ 50meV. 



ity CTint/int I On the other hand, is almost unaffected by 
the VC, retaining the characteristic maximum at about 
2 A (labeled as 2 A maximum in the following). The in- 
dependencies of the energies of the B. A. mode and of the 
2A-peak in the local conductivity are shown in Fig.[T2ja). 
The energy of the B. A. mode can be seen to be propor- 
tional to t±, except for the region of high values of t±] 
that of the 2A-peak is approximately tj^-independent. 

The 2A-maximum appears also in the total c-axis con- 
ductivity, as shown in Fig. [TlTb). The spectral weight 
of the corresponding peak (labeled as T2 in the follow- 
ing) is proportional to t'^. As documented in Fig. [T^ b). 
for low values of tx, Ti dominates and T2 cannot be re- 
solved. For high values of t±, on the other hand, T2 is 
the most pronounced feature since the B. A. mode merges 
with the continuum and Ti can hardly be resolved. Both 
features can be seen for intermediate values of t±, e.g., 
t± = 70meV. The ^x-dependencies of the energies of the 
structures Ti and T2 in the total c-axis conductivity are 
given in Fig. [TWa). 



C. Comparison with experiment and with the 
theory proposed by Shah and Millis 

The c-axis conductivity displays two 
superconductivity- induced structures (modes): Pi 
and P2 in the experimental data, and Ti and T2 in the 
theoretical spectra. In what follows, we argue that the 
features Pi and P2 can be attributed to Ti and T2, 
respectively. 

First, we summarize the relevant trends of the struc- 
tures Pi and P2 as observed in the experimental data of 
bilayer compounds, in particular Y-123 and related sys- 
tems. Some of the trends are demonstrated in Fig. [131 
(El) The frequency of Pi increases with increasing hole 
concentration pi^^Oiis^iSi^^ 

(E2) The spectral weight (SW) of Pi first increases with 
increasing p, then saturates for p w 0.12, and for higher 
values of p, Pi broadens and its SW gradually decreases; 
for p > 0.15, Pi cannot be resolved anymore. ^-2, 
(E3) The frequency of Pi decreases when going - for a 
fixed doping level - from Y-123 over Nd-123 to La-123.-^ 
In this sequence of materials, the distance between the 
closely-spaced planes increases. 

(E4) In the YPr-123 system, the doping level can be mod- 
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FIG. 13; (a) Doping dependence of the diflterence Keac{T <C 
Tc) ~ Reac{T ^ Tc) for Y-123. The abbreviations UD, OPT, 
and OD stand for underdoped, optimum doped, and over- 
doped. The values of p are 0.093, 0.116, 0.124, 0.155, and 
0.194. Also shown are the data for Yo.86Cao.i4-123 with 
Tc — 84 K and p — 0.176. Details concerning the samples 
and the experiment are given in Refs. [20l |37| . (b) The spectra 
of the difference for the sequence R — 123 {R = Y, Nd, and 
La) with p ~ 0.12. In this sequence, the distance between the 
Cu02 planes within a bilayer increases. The inset shows the 
original spectra for Nd-123. Adapted from Fig. 2 of Ref . [20l. 



ified either by changing the oxygen concentration or by 
partially replacing Y with Pr. The replacement leads to 
a decrease of p. By combining the two approaches, it is 
possible to obtain various combinations of p and the NS 
dc-conductivity along the c-axis, a^c- For a fixed doping, 
there is no pronounced correlation between the SW of Pi 
and adc^ 

(E5) The structure P2 can be resolved only for p > 0. 10*2^ 
(E6) The frequency of P2 slowly decreases with increas- 
mg p. 

(E7) The SW of P2 increases with increasing p^ 
(E8) In the YCa-123 system, it is possible, similarly as 
in the case of YPr-123, to obtain various combinations 
of p and (Tdc- For a fixed doping, the SW of P2 increases 
with increasing o'dc'^ 

The distance between the closely-spaced planes can be 
expected to be correlated with the strength of the intra- 
bilayer electronic coupling. The observation (E3) thus 
suggests a relation between Pi and the coupling. Fur- 
ther, the dc-conductivity is likely to reflect the strength 
of the coupling through the spacing layer. The observa- 
tions (E4) and (E8) thus seem to imply an independence 



of Pi and a dependence of P2 on this coupling. 

Second, we review the corresponding properties of Ti 
and T2 resulting from our computations. The labels 
(T1)-(T8) parallel those used in the summary of experi- 
mental facts. 

(Tl) The energy of Ti increases with increasing t±, see 
Figs.iniandlini 

(T2) The SW of Ti first increases with increasing t±, 
then saturates for t± w 50 meV, and for higher values of 
t±, Tl broadens and gradually vanishes, see Fig. [TUf f). 
The broadening is due to the fact that Ti reaches the 
continuum. 

(T4) For a given t±, Ti does not change significantly with 
increasing t'j_ (not shown). Note that t'j_ determines adc- 
(Tdc is approximately proportional to t'^. 
(T6) The energy of T2 is approximately equal to 2A for 
low values of t± and somewhat larger than 2 A for higher 
values of t± of the order of 100 meV, see Fig. [T^ aV 
(T7) The SW of T2 increases with increasing t± (not 
shown) . 

(T8) The SW of T2 also increases with increasing t'^ , see 
Fig. HH For a given t± , the SW of T2 is approximately 
proportional to i'^, similarly as (Tdc. 

A comparison between the items (E1)-(E4) and (Tl)- 
(T4) suggests that Pi could be attributed to Ti provided 
that tj_ increases with increasing p. This crucial assump- 
tion is consistent with results obtained using the bilayer 
t~J model and the Gutzwiller approximation42ii^ Alter- 
natively, it can also be justified considering a pseudogap 
(PG) competing with superconductivity and the reported 
p- and fe-dependencies of the magnitude of the PG and of 
the coherence peaks due to superconductivity4i»^ With 
decreasing p, the magnitude of the PG increases and the 
area of the part of the Brillouin zone dominated by the 
PG, centered around the antinode, expands. At the same 
time, the area of the part with pronounced Bogolyubov 
quasiparticles shrinks towards the BZ diagonal. The im- 
portant point is that at the BZ diagonal, t± is probably 
the smallest)^ The shrinkage of the area of strong super- 
conducting correlations might thus lead to a decrease of 
an effective t±_ determining the energy of the low-energy 
mode. 

The properties (E5)-(E8) are in agreement with at- 
tributing P2 to T2: The items (E5) and (E7) can be un- 
derstood in terms of (T7) and (T8), when combined with 
the assumption of t± ^ p and with the obvious fact that 
t'j_ increases with increasing oxygen concentration; (E6) 
is consistent with (T6), when combined with the exper- 
imental fact that, around optimum doping, A decreases 
with increasing p. Further, (E8) can be understood based 
on (T8). 

An obvious problem of the proposed assignment of Ti 
to Pi and T2 to P2 is that t^(Ti) > lu{T2) for the rele- 
vant range of parameters (see Fig. I12p . whereas for un- 
derdoped materials with 0.10 < p < 0.15 the frequency 
of Pi is lower than that of P2i^ The reason for this 
discrepancy is probably the following: Experimentally, 
Ci;(P2) « 100 meV which requires A w 50 meV, and such 
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FIG. 14: Figure illustrating the potential of the model to 
provide the order of the spectral structures consistent with 
experimental data on underdoped Y-123. (a) Real part of the 
local conductivity (Jbi/bi consisting of the result of the BCS- 
approach with A = 50meV and a broad Lorentzian described 
in the text, (b) Real part of the total c-axis conductivity 
obtained using the local conductivity (Tbi/bi shown in (a) and 
the other three local conductivities resulting from the BCS- 
approach. A constant value of the ratio t'xmax/^Xmax of 0.2 
has been used. 



a high value cannot be achieved using the present self- 
consistent theory with reasonable values of input param- 
eters. The effective gap may have a contribution due to 
the pseudogap, that is not included in the present theory. 
Within the BCS approach, A is an input parameter. Un- 
fortunately, the resulting spectra of cjc for the interesting 
region of the parameter space (A > t±) contain over- 
lapping resonances that cannot be easily disentangled. 
The complication is due to the absence of an incoher- 
ent background in the local conductivities. Motivated by 
this observation, we have supplemented the component 
fbi/bi with a broad Lorentzian, —iA/{ujj^ — — io^r^), 
UJL = 0.4 eV, Tl = 0.15 eV, and A ~ t'^. The t±- 
dependence of the spectral-weight parameter A is consis- 
tent with Eq. (PT|) . The results are shown in Fig. [TD It 
can be seen that the i^-dependence of the total conduc- 
tivity shown in (b) resembles the doping dependence of 
the data^", including the interplay of Ti and T2: For low 
values of t± , the spectra are dominated by Ti ; for inter- 
mediate values of 20 — 30 meV, both spectral structures 
are present, and with increasing t±, Ti gradually hides 
in the continuum. 

The discrepancy in the order of the spectral structures 
does not occur for optimally doped and overdoped Y-123 
samples, where Pi is buried in the continuum part of the 
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FIG. 15: Representative examples of the spectra of Re CTc com- 
pared with the conductivities of a related single- layer system, 
and with the spectra of the same quantity obtained using 
the approximation proposed by Shah and Millis. (a) Real 
part of the total c-axis conductivity calculated using the spin- 
fermion model with txmax ~ 250 meV, tYmax = 75 meV, and 
with txmax = 30 meV, t'xmax ~ meV. The solid and the 
dashed lines correspond to the superconducting and the nor- 
mal state, respectively. The spectra for txmax = 250 meV and 



— 75 meV on an extended scale are shown in the inset. 



(b) Real parts of the in-plane conductivity and of the c-axis 
conductivity for a single-layer superconductor with the plane 
spacing of dbi + rfint and with the c-axis hopping parameter of 
the form of Eq. (0 and the maximum of 75meV. (c) The same 
as in (a) but using the approach proposed by Shah and Millis, 
where the nondiagonal components of the conductivity (Tbi/int 
and (Tint/bi are neglected, and the diagonal component CTbi/bi 
(o"int/int) IS approximated by the conductivity of the single- 
layer superconductor with the hopping parameter equal to t± 
{t'j_). Only the superconducting state spectra are shown. 



spectra located above the maximum P2 . The experimen- 
tal spectra (see e.g. Fig. 1 of Ref. IT3 ) are similar to the 
calculated ones corresponding to the values of t± of a few 
hundreds of meV. As an example, we show in Fig. llSf al 
the NS and SC-state spectra of Reo-c corresponding to 
t±ma.K = 250 meV and t'^^^^^ = O.Si^max- The maximum 
T2 at about 70meV (600 cm^^) probably corresponds to 
P2 occurring at a slightly higher frequency in the exper- 
imental spectra. 

Next we compare our calculated spectra of Re CTc with 
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the conductivities of a related single-layer superconduc- 
tor and with those computed along the lines of the the- 
ory proposed by Shah and Millis (SM) .— Figure fTSFa) 
shows, besides the spectra of Rbctc (ij_max = 250 meV, 
t±ma.x — 75 meV) discussed above, those corresponding 
to tj_max — 30 meV, i±max = mcV. The former (the 
latter) represent the case of strong (weak) intra-bilayer 
coupling. The corresponding renormalized values of the 
(normal-state) bilayer splitting are ca 80 and 10 meV. 
Note that the former value is close to that of Bi-2212 
as obtained by photoemission experiment o^^'^^i"'^^'^^ . As 
discussed in the context of Figs. [TT] and [T^] only one peak 
is present in the superconducting state spectra: T2 in 
the former case and Ti in the latter. It is instructive to 
compare the SC-state conductivities with the solid line 
of Fig. I15r b) representing the c-axis conductivity of a 
model single-layer superconductor described in the cap- 
tion. This allows one to identify the features specific to 
the bilayer compounds: (i) the peak Ti (for small val- 
ues of t±), (ii) the peak T2 (for large values of t± and 
t'j_ 0), and (Hi) a hump in the mid-infrared (for large 
values of t±). For t^max = 250 meV, the hump is cen- 
tered at 500 meV. It has the same origin as the 160 meV 
maximum in Fig. [3] and it can be attributed to the (up- 
per) plasma mode of the bilayer unit, see the discussion 
following Eq. ([50)) . Both Ti and T2 would appear for 
intermediate values of t±, as discussed in the context of 
Fig. [12] and Fig. [141 The broad band centered around 
120 meV in Fig.fTSTa') appears also in the conductivity of 
the single-layer superconductor and is thus not specific 
to the bilayer compounds. 

Also shown in Fig. I15f b) is the in-plane conductivity 
of the single-layer superconductor. It can be seen that 
the shapes of the two conductivities are fairly similar. 
The main differences are: (a) for T = 20 K, a narrow 
Drude term can be resolved only in Reuab, and (b) the 
onset of RecTc is more gradual than that of Keaab- Both 
(a) and (b) are due to the difference between the matrix 
element of aab, *-e-, the in-plane quasiparticle velocity, 
and t±. Both aat and CTc exhibit a maximum around 
110 meV and both decrease with increasing energy above 
this maximum. The origin of the maximum has been 
addressed in Refs. [s^, H^. The calculated c-axis conduc- 
tivity of the single-layer superconductor is qualitatively 
similar to the measured conductivity of La2-2;Srj;Cu04 
reported in Ref. H^. 

Figure [Tbtc) shows the c-axis conductivities of the 
same bilayer systems as in (a), but calculated using the 
approach proposed by SM, where the nondiagonal com- 
ponents of the conductivity tJbi/int and dint/bi are ne- 
glected, and the diagonal component <Jt,\/h\ (cint/int) is 
approximated by the conductivity of the single-layer su- 
perconductor with the hopping parameter equal to t±_ 
[t'^ ) . When calculating the latter conductivities, SM fur- 
ther replace the fcz-dependent Green's function with a 
fc^-independent one of a two-dimensional model. Instead 
of using this approximation, we have obtained the con- 
ductivities of the model single-layer superconductors by 



scaling the one shown in Fig. [TSl' b). We have checked 
that for the present values of the parameters, the results 
of the two approaches are almost the same. By compar- 
ing panels (a) and (c) of Fig.[Tni we easily identify differ- 
ences between the SM theory that does not involve the 
bilayer splitting and our improved approach, that does. 
For small values of t±_ , the SM theory provides a peak in- 
side the gap, corresponding to a bilayer plasmon, similar 
to the transverse Josephson plasmon of the phenomeno- 
logical modelii^ The Ti peak of our theory is located at a 
slightly higher energy and its interpretation is different. 
What it has in common with the transverse Josephson 
plasmon is that both are associated with oscillations of 
the relative phase of the two closely-spaced planes. For 
high values of t± , our approach yields the pair-breaking 
peak T2 absent at the SM level. At high energies (above 
100 meV), the results of the two approaches are similar. 

Interestingly, the NS spectra corresponding to t_Lmax = 
250 meV shown in Fig.fTSTa) do not display any clear sig- 
nature of the bilayer splitting, consistent with experimen- 
tal data. The contribution of the bonding-antibonding 
transitions is hidden in the mid-infrared region. This is 
because the renormalization of the quasiparticles leads to 
a very broad absorption band [see Fig. [TUTa)] shifted to- 
wards the mid-infrared both by the the vertex corrections 
and by the Coulomb effects [see Fig. [TUfb) and (e)]. 

Experiments reveal an increase of the optical spectral 
weight in the far-infrared below Tc, with a possible in- 
terpretation in terms of a decrease of the c-axis kinetic 
energy associated with the superconducting transition.— 
The picture resulting from our calculations, restricted to 
the case of small values of i_Lmax and t'^^^^ = 0, is the 
following: Below Tc, the peak Ti forms, gaining spectral 
weight from a broad interval of energies, and the spec- 
tral weight at low energies increases. The total spectral 
weight, proportional to the negatively taken effective ki- 
netic energy K\^i, decreases, provided that xsf (below 
Tc)— XSF (above Tc). However, changes of xsf upon 
entering the superconducting state, in particular the for- 
mation of the resonance mode, can lead to a slight in- 
crease of the total spectral weight and the corresponding 
decrease of K\^\, the mechanism being connected to that 
outlined in Ref. The issue is fairly complex and will 
be addressed in a separate publication. 



IV. SUMMARY AND CONCLUSIONS 

We have constructed a realistic microscopic model of 
the c-axis infrared response of bilayer cuprate supercon- 
ductors allowing us to interpret the superconductivity- 
induced modes occurring in the experimental data. 

For the simpler case of insulating spacing layers, the 
local conductivity of the intrabilayer region does not pos- 
sess a condensate contribution in the superconducting 
state [S{uj) in Recr] as assumed within the phenomeno- 
logical Josephson superlattice model. Instead, it displays 
a collective mode at a finite frequency that is proportional 
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to the interplane hopping amphtude t±. This has been 
shown to be a consequence of the gauge invariance. The 
nature of the mode is similar to that of the Bogolyubov- 
Anderson mode that participates in the longitudinal re- 
sponse of a homogeneous superconductor. It is associated 
with charge oscillations between the planes and, for small 
values of t±, also with oscillations of the relative phase 
of the two planes. This physical picture is fairly simi- 
lar to that of the transverse plasmon of the Josephson 
superlattice model. In the total c-axis conductivity the 
mode is shifted towards higher energies by the interplane 
Coulomb interaction. 

A nonzero amplitude of the hopping through the spac- 
ing layer implies a finite conductivity of this layer. This 
local conductivity exhibits a peak at a frequency slightly 
higher than 2Aniax, that can be interpreted as a pair 
breaking peak. The simple picture is such that two Bo- 
golyubov quasiparticles are involved: one from the bond- 
ing band and the other from the antibonding. The rea- 
son, why the peak appears in the inter-bilayer conduc- 
tivity and not in the intra-bilayer one, is the following. 
Excited states behind the former (latter) conductivity 
involve quasiparticles from different bilayers (from the 
same bilayer). Only the latter are thus strongly modi- 
fied by including the final-state interactions, that are re- 
stricted to individual bilayers. The peak permeates into 
the total c-axis conductivity. 

A series of arguments has been presented assigning 
the collective mode to the low energy superconductivity- 
induced mode of the experimental data (interpreted pre- 
viously in terms of the Josephson superlattice model) and 
the pair-breaking maximum to the superconductivity- 
induced peak centered around 1000cm~^. The arguments 
concern the doping dependence of the frequencies and the 
spectral weights of the peaks and the impact of various 
substitutions. A crucial assumption, connecting the the- 
ory and the experiment, is that the effective t± decreases 
with decreasing doping. The trends of the underdoped 
regime, in particular, the appearance of the collective 
mode below Tc, the increase of its frequency with increas- 
ing doping and its disappearance below optimum doping 
can all be reasonably reproduced and understood using 
this assumption. Admittedly, the values of tj_ of a few 
tens of meV needed to fit the data are smaller than those 
deduced from photoemission experiments. The main fea- 
tures of the data of optimally doped Y-123 can be rea- 
sonably reproduced with i j.max of 250 meV, which corre- 
sponds to the maximum distance between the renormal- 
ized bands of ca 80 meV. 
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APPENDIX A: BCS LEVEL OF THE THEORY 

In this appendix, we show several results obtained at 
the BCS level, where extensive analytical simplifications 
can be made. For the sake of brevity, we restrict ourselves 
to the case of t^^^^ = 0. We employ the BCS interaction 
of d-wave symmetry which acts in the individual Cu02 
planes (labeled as 1 and 2) within a bilayer unit. The 
corresponding Hamiltonian reads 

^^BCS = ^BkC^Bka^^Bka + X! ^ ^f'^^Aka'^ Ak,y + 

ka kcr 

+ X! 2Vfcfc'cJj^|cJj _j,j^c,j _^,j^c^^,^ , (Al) 

fcfc',riG{l,2} 

where Vkk' is introduced in the main text. The factor of 
2 is for later convenience. The interaction term can be 
written as 

^ Vkk'iBBBB + BBAA + AABB + AAAA+ 



kk' 



BAB A + BAAB + ABAB + ABBA) . (A2) 



where BBBB, e.g., stands for c^^^c^ -ki'^B -k'l^Bk'] - 
The first four terms in Eq. (|A2[) provide the pairing in- 
teraction, the other four terms play an important role 
in the vertex corrections. The above pairing interaction 
is equally distributed among the two symmetry chan- 
nels and produces the same gap in both bands, hence 
— AfcTi with Afc determined by Eq. pT|) . 
Because the selfenergy depends on k only, we can sum 
over the Matsubara frequencies explicitly. The evaluation 
of the NV response function given by Eq. (|2ip leads to 



^A/B(k 



ttNV(i), n t ^ e dbi ^ 



±k' 



li[l~nF{EAk)-nF{EBk)] 



+l2[nF{EAk)~nFiEBk)] 



hiv + E^+iS huj-E^ + iS, 



fiuj + Ef^ + i5 huj — Ep, + i6 ^ 



where Ej^ — EAk^Esk and I1/2 are the coherence factors 

SAk^Bk + ^Ak^Bk 



h/2 



EAk Esk 



(A4) 



Note that in the limit of t±f, 0, the factors li as well 
as nF{EAk) — nF^EBk) vanish and Ilbi-bi becomes zero. 



(A3) 
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For the same reason, there is no intraband contribution 
dm) in the BCS case. 

The renormahzed vertices Tab and Tba satisfy the 
Bcthe-Salpetcr equations similar to Eq. now con- 

taining the BCS interaction: 



TABik, iE, ihv) = t_Lfcro 



kBT 
N 



k' ,iE' 

T3 Geik', lE') TAB{k\ iE' , ihv)gA{k' , iE' + ihvY 



GAik', iE') TBA{k', iE', ifiv)gB{k' , iE' + ihv) 



T3 . 

(A5) 



Since the usual BCS interaction couples with T3 , not with 
To like the spin-fluctuations, the above equation contains 
the additional matrices t^. The corresponding equation 
for the vertex Tba differs from Eq. (jASP in the sign 
of the t±kTQ term only. By inserting the form of V^fc' 
explicitly, one finds, that the vertices can be cast to: 
T AB{k,ifiv) = t±kTo + \wkC{ihv) and T BA{k,ifiv) — 
—t±kTo -\- XwkC {ihi>) . The quantity C{ihv) is a 2 x 2 ma- 
trix, independent of iE because of the non-retarded BCS 
interaction and independent of k because of the separa- 
ble form of Vfcfc' . We express it as a linear combination 
of the Pauli matrices: C{ihv) — J2a=o '''aC a {itiv)- The 
Bethe-Salpeter equation (|A5[) can be then converted to 
a linear system of equations for the coefficients Ca{ihv) 

VaCa{ihv)-^ Aaf){ifLiy) Cpiifiv) = <^aa{iTT-v) , a,P = 0, 

(A6) 



with 77c, = -1-1 (a = 0,3), r?a = -1 (a = 1,2). The 
coefficients of the system are given by 



^Tq$q;3 = '^"KT 51 Wkt^k[Q B{k,iE) Tp g A{k,iE+ifiv) 



k.iE 



gAik,iE)T(3gBik,tE + ihiy)] (AT) 



and 



TqAo^ = A—— ^ wi [gB{k, iE) Tp gA{k, iE+ifiv) 



N 



kAE 



+ gA{k,iE)TpgB{k,iE + ihiy)] . (A8) 

Finally, to get the response function, wc insert the renor- 
mahzed vertices into (|23p and with the help of Trr^ = 
25ao obtain 



bi-bi - ^^bi-bi - -TT-T y2 '^0/30/3 . (A9) 



ttVC _ ttNV 

r^Ki—Ki — ri 



All the above equations can be analytically continued to 
the real axis explicitly. For each frequency required, we 
have to evaluate the coefficients and Ac,^, solve the 
linear system (|A6P and find the response function (jA9p . 
Thanks to the simple form of the selfenergy, the Matsub- 
l,2^a summations can be again performed analytically, as 
in Eq. (jA3p . but lead to more cumbersome expressions 
due to the additional Tp matrix in Eqs. IA7l andf 
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